This section is intended to introduce various aspects of the art, which may be associated with aspects of the disclosed techniques and methodologies. A list of references is provided at the end of this section and may be referred to hereinafter. This discussion, including the references, is believed to assist in providing a framework to facilitate a better understanding of particular aspects of the disclosure. Accordingly, this section should be read in this light and not necessarily as admissions of prior art.
Reservoir simulators solve systems of equations describing the flow of oil, gas and water within subterranean formations. In a reservoir simulation model, the subterranean formation is mapped on to a two- or three-dimensional grid comprising a plurality of cells. Each cell has an associated equation set that describes the flow into the cell, the flow out of the cell, and the accumulation within the cell. For example, if the reservoir is divided into 1000 cells, there will be 1000 equation sets that need to be solved.
To model the time-varying nature of fluid flow in a hydrocarbon reservoir, the solution of the equations to be solved on the grid cells varies over time. In the reservoir simulator, solutions are determined at discrete times. The time interval between solutions is called the timestep. For example, the reservoir simulator may calculate the pressures and saturations occurring at the end of a month, so the timestep is a month and a single solution to the equation set is needed. To calculate the changes in pressure and saturation over a year, the simulator in this example calculates twelve monthly solutions. The time spent solving this problem is roughly twelve times the time spent solving the single month problem.
The size of the timestep that a simulator can take depends on a number of factors. One factor is the numerical method employed to find the solution. As some reservoir models may have thousands or millions of cells, various methods have been proposed to efficiently solve the large number of equations to be solved by a reservoir simulation model. One known strategy for finding the solutions to these systems of equations is to use an iterative root-finding method. These methods find approximate solutions that get progressively closer to the true solution through iteration and solution updating. Newton's Method is one type of iterative root-finding method in general use. In Newton's Method, the set of simulation equations are cast into a form that makes the solution an exercise in finding the zeros of a function, i.e. finding x such that ƒ(x)=0.
FIG. 1 is a graph 8 that depicts Newton's Method for a single equation. Curve 10 is the function ƒ(x). What is sought is the value x where ƒ(x)=0, indicated by point 12. The initial guess is x0. The second guess is calculated by taking the line 14 tangent to ƒ(x) at x0 and applying the formula x1=x0−(ƒ(x0)/ƒ′(x0)). Here ƒ′(x) denotes the derivative of the function ƒ(x) and is the slope of the tangent line at x. The third guess, x2, uses the line 16 tangent at the second guess (x1) and applies the same formula, x2=x1−(ƒ(x1)/ƒ′(x1)). Continuing this iterative algorithm one can get very close to the root of ƒ(x), i.e., point 12, in a modest number of iterations. If the curve ƒ(x) is well-behaved, quadratic convergence can be achieved.
Reservoir simulators have expanded Newton's Method to solve for the many thousands of equations at each timestep. Instead of one equation a system of equations is used:
                                                        f              1                        ⁡                          (                                                x                  1                                ,                …                ⁢                                                                  ,                                  x                  n                                            )                                =          0                ⁢                                  ⁢                                            f              2                        ⁡                          (                                                x                  1                                ,                …                ⁢                                                                  ,                                  x                  n                                            )                                =          0                ⁢                                  ⁢        ⋮        ⁢                                  ⁢                                            f              n                        ⁡                          (                                                x                  1                                ,                …                ⁢                                                                  ,                                  x                  n                                            )                                =          0                                    [                  Equation          ⁢                                          ⁢          1                ]                where ƒ1 (x1, . . . , xn)=0 are the reservoir simulation equations for grid block 1 containing the variables x1 through xn, and n is the number of grid cells. The variables, x1, . . . , xn, are typically pressures and saturations at each cell.
To apply Newton's method to this system of equations the tangent of the function is needed at each iteration, like that described above for the single equation above. The tangent of this matrix A is called the Jacobian J and is composed of the derivatives of the functions with respect to the unknowns.
                    J        =                  [                                                                                          ∂                                          f                      1                                                                            ∂                                          x                      1                                                                                                  …                                                                                  ∂                                          f                      1                                                                            ∂                                          x                      n                                                                                                                          ⋮                                            ⋱                                            ⋮                                                                                                          ∂                                          f                      n                                                                            ∂                                          x                      1                                                                                                  …                                                                                  ∂                                          f                      n                                                                            ∂                                          x                      n                                                                                                    ]                                    [                  Equation          ⁢                                          ⁢          2                ]            
As in the case with one equation an initial guess {right arrow over (x)}0 is made, where {right arrow over (x)}0 is a vector of solutions. Each subsequent guess is formed in the same manner as that for a single variable, where {right arrow over (x)}i is formed from the previous guess {right arrow over (x)}i-1 using the following equation:{right arrow over (x)}i={right arrow over (x)}i-1−(f({right arrow over (x)}i-1)/J({right arrow over (x)}i-1))  [Equation 3]This equation can be rewritten asJ({right arrow over (x)}i-1)({right arrow over (x)}i−{right arrow over (x)}i-1)=−f({right arrow over (x)}i-1)  [Equation 4]which is a matrix equation of the form A{right arrow over (x)}={right arrow over (b)}. The solution is thought to be converged when either the term ({right arrow over (x)}i−{right arrow over (x)}i-1) or −f({right arrow over (x)}i) approaches zero, i.e. is below a small threshold, epsilon (ε). This idea is applied to each of the thousands or millions of cells over hundreds of timesteps in a reservoir simulator.
FIG. 2 is a flowchart 20 showing the steps of Newton's Method for a system of equations. At block 21, a solution vector {right arrow over (x)}i, representing the solutions for the system of equations, is set to an initial guess {right arrow over (x)}0. At block 22 a Jacobian Matrix Ji and a vector {right arrow over (b)}i is constructed for all cells Zn, or equation sets associated with the cells, using the solution vector {right arrow over (x)}i. At block 23 a new solution estimate {right arrow over (x)}i1 is obtained for all cells Zn. At block 24 it is determined how many cells are not converged, which may be defined as having associated equation sets that have not satisfied a convergence criterion. If the number of unconverging cells is zero, then at block 25 the method stops. However, if at block 26 it is determined that the number of unconverging cells is not zero, then at block 27 the solution guess vector {right arrow over (x)}i is set to the new solution estimate {right arrow over (x)}i1, and the process returns to block 22. The method repeats until all cells are converged. At that point the system of equations may be considered solved.
While the Newton's Method provides a simple way to iteratively solve for solutions to systems of nonlinear equations, its effectiveness is lessened when the equations in a system of equations do not uniformly converge. For example, a two-dimensional grid 30 is shown in FIGS. 3A, 3B, 3C, and 3D as having 169 cells. Each cell has an equation or equation set associated therewith. Prior to any iteration (FIG. 3A) the 169 equations are used as inputs to Newton's Method. After one iteration of Newton's Method (FIG. 3B) the equations related to 111 cells have converged (shown by the lighter-colored cells 32). In other words, some areas of the reservoir have found solutions such that the term ({right arrow over (x)}i−{right arrow over (x)}i-1) is below ε for those areas. The unconverged cells are shown as darker colored cells 34. The second iteration (FIG. 3C) uses all 169 equations, and the equations related to 147 cells have converged thereafter (as indicated by reference number 36). After the third iteration (FIG. 3D)—which also uses all 169 equations—is completed, all equations have converged. This example highlights a drawback of Newton's Method: even though most equations have converged to a solution after a single iteration, Newton's Method uses the equations for all cells for all iterations. In other words, the size of the system of equations to be solved does not change after each iteration. For a large two- or three-dimensional reservoir simulation grid having thousands or millions of cells, the time and computational power required to perform Newton's Method repeatedly may be prohibitive.
In general, the linear system constructed by Newton's method is very large and expensive to solve for typical reservoir simulation models with complex geological structures and difficult physical properties. While maintaining numerical stability, fully implicit schemes employing the Newton's Method require substantial CPU time and a large memory footprint. More explicit schemes are cheaper in CPU time per timestep, but have difficulty taking viable time steps with reasonable sizes, given their stability limits. The Adaptive Implicit Method (AIM) is a natural choice to balance implicitness/stability and CPU/memory footprint. The basic concept of AIM is to combine multiple formulations with different degrees of implicitness. Therefore, the simulation can use a fully implicit formulation (e.g., solving pressure and saturations at the same time) to maintain unconditional stability for cells with constraining stability, while using cheaper and more explicit formulations (e.g., solving pressure only, explicitly updating saturation solutions afterwards) for the remainder of the reservoir. This type of strategy is known as IMplicit Pressure Explicit Saturation (IMPES). AIM reduces the size of the Jacobian matrix without sacrificing the numerical stability and timestep sizes. AIM relies heavily on robust stability estimation and prediction to determine how to distribute reservoir simulation cells into different formulations, i.e., potentially unstable cells in the implicit region, and less stable cells in the explicit region.
FIGS. 4A-4C show an example of how the Adaptive Implicit Method may be applied in a reservoir simulator. A 17×17 reservoir simulation grid 40 is depicted at a particular timestep, demonstrating a five-spot water injection pattern with a water injector at the upper-left corner 41 and a producer at the lower-right corner 42. FIG. 4A is the map of water saturation at the given timestep with arrow 43 indicating the direction of the water front movement. The degree of shading of each cell indicates the amount of water saturation therein. FIG. 4B graphically illustrates relative stability at each cell at the beginning of the timestep with shaded cells 44 denoted as unstable. These cells 44 are assigned to be solved implicitly. The stability of each cell will be recalculated at the end of the timestep, which is shown in FIG. 4C. As the water front moves during the timestep, the stability at each cell may change, and this is shown in FIG. 4C as a different group of cells 45 being denoted as unstable. This group of unstable cells 45 is assigned to be solved implicitly for the next timestep.
A typical AIM scheme will preset the implicit/explicit formulations at the beginning of each timestep using the stability information calculated from previous timestep. Although the reservoir model contains multiple formulations, the formulation on any given cell will remain the same throughout the timestep, This approach has a few potential issues. First, all the stability criteria have been derived from linearized systems. For highly nonlinear systems, the stability criteria are not very robust and reliable, e.g., as in thermal recovery processes. Second, the stability limits are calculated at the beginning of each timestep. The situation often arises that a cell was stable according to the stability calculation and put into the explicit region, but by the end of the timestep the cell is not stable any more because of the nonlinear nature of the problem, as shown in FIGS. 4B-4C. For example, one of the commonly used stability criteria is expressed as the CFL (Courant-Friedrichs-Lewy) number, and the stability condition for IMPES formulation is CFL<1. In a simplified two-phase oil-water simulation model, the CFL number can be expressed as
                              C          ⁢                                          ⁢          F          ⁢                                          ⁢          L                =                                                            q                T                            ⁢              Δ              ⁢                                                          ⁢              t                                      V              p                                ⁢                                    ⅆ                                                F                  w                                ⁡                                  (                                      S                    w                                    )                                                                    ⅆ                              S                w                                                                        [                  Equation          ⁢                                          ⁢          5                ]            where Δt is the timestep size, qTΔt/Vp the volumetric throughput, and Fw(Sw) is the fractional flow, which is a function of water saturation Sw according to
                                          F            w                    ⁡                      (                          S              w                        )                          =                              λ            w                                              λ              w                        +                          λ              o                                                          [                  Equation          ⁢                                          ⁢          6                ]            where λw, λo are water and oil mobility respectively. If a cell has a low CFL number (less than 1, for example) at a given timestep, it is assumed that cell is stable, and an explicit formulation is used to solve the flow equations at that timestep. If a cell has a high CFL number (greater than 1, for example) at a given timestep, it is assumed the cell is unstable, and an implicit formulation is used to solve the flow equations at the timestep. FIG. 5 shows a graphic representation of how the CFL number is calculated for a simulation cell. Fw′(Sw) is the derivative of the function Fw(Sw) 50 which is the slope of the tangent line at any given point. With a given timestep size and fixed volumetric throughput, the CFL number is proportional to this derivative as expressed in Equation 5. The fractional flow Fw1(Sw1) at the beginning of the timestep (with water saturation Sw1) has a CFL number just within the stability bound, represented by the relatively flat slope of tangent line 51, and dictates an explicit formulation. The fractional flow Fw2 (Sw2) at the end of the timestep (with water saturation Sw2) has a CFL number much larger than the stability bound, represented by the steep slope of tangent line 52, and may cause a numerical stability issue.
Several strategies may be used to remedy this stability issue, but each strategy imposes unwanted costs on the simulation. A conservative CFL limit may be used, but this may impose a more stringent timestep constraint on the simulation run, or distribute too many simulation cells into the implicit formulation. The timestep may be rerun with the unstable cells being switched to the implicit formulation, but for large and highly nonlinear systems this approach will inevitably incur unnecessary calculations and slow down the simulation. A timestep cut could be triggered to satisfy the stability limit, but this would increase the number of calculations to solve the simulation. Finally, a small amount of unstable cells may be permitted to be run in the simulation, assuming that the local instability will dissipate and disappear eventually. This approach sacrifices numerical stability for simulation runtime performance.
Another potential inefficiency of the AIM method could be aggravated if a conservative CFL limit is used. A cell might require a fully implicit scheme at the beginning of a timestep. However, the particular cell could quickly converge to a solution that is well within the stability limit. As FIG. 5 shows, the CFL number at the end of another timestep (with water saturation Sw3), represented by the relatively flat slope of tangent line 53 could be much lower than the beginning of the timestep (with water saturation Sw2), represented by the steep slope of tangent line 52. Assigning the associated simulation cell to an implicit formulation scheme will likely over-constrain either the timestep size or the computational scheme.
Various attempts have been made to reduce the time required to solve a system of equations using implicit or adaptive implicit methods. Examples of these attempts are found in the following: U.S. Pat. Nos. 6,662,146, 6,052,520, 7,526,418, 7,286,972, and 6,230,101; U.S. Patent Application No. 2009/0055141; “Krylov-Secant Methods for Accelerating the Solution of Fully Implicit Formulations” (SPE Journal, Paper No. 00092863); “Adaptively-Localized-Continuation-Newton: Reservoir Simulation Nonlinear Solvers that Converge All the Time” (SPE Journal, Paper No. 119147); “Preconditioned Newton Methods for Fully Coupled Reservoir and Surface Facility Models” (SPE Journal, Paper No. 00049001); Gai, et al., “A Timestepping Scheme for Coupled Reservoir Flow and Geomechanics on Nonmatching Grids” (SPE Journal, Paper No. 97054); Lu, et al., “Iteratively Coupled Reservoir Simulation for Multiphase Flow” (SPE Journal, Paper No. 110114); and Mosqueda, et al., “Combined Implicit or Explicit Integration Steps for Hybrid Simulation”, Earthquake Engineering and Structural Dynamics, vol. 36 no. 15, pp 2325-2343 (2007). While each of these proposed methods may reduce the time necessary to solve a system of equations, none of the methods reduce the number of equations (or cells) required to be solved using an implicit method employing Newton's Method. Furthermore, none of the references disclose a change in formulation method during a timestep. What is needed is a way to reduce the number of equations needed to be solved in successive iterations of an iterative solver.